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Abstract 

Here we use a very simple conceptual model in an attempt to reduce essential parts of 
the complex nonlinearity of abrupt glacial climate changes (the so-called Dansgaard- 
Oeschger events) to a few simple principles, namely (i) the existence of two different 
climate states, (ii) a threshold process and (iii) an overshooting in the stability of the 
system at the start and the end of the events, which is followed by a millennial-scale 
relaxation. By comparison with a so-called Earth system model of intermediate com- 
plexity (CLIMBER-2), in which the events represent oscillations between two climate 
states corresponding to two fundamentally different modes of deep-water formation in 
the North Atlantic, we demonstrate that the conceptual model captures fundamental as- 
pects of the nonlinearity of the events in that model. We use the conceptual model in or- 
der to reproduce and reanalyse nonlinear resonance mechanisms that were already sug- 
gested in order to explain the characteristic time scale of Dansgaard-Oeschger events. In 
doing so we identify a new form of stochastic resonance (i.e. an overshooting stochastic 
resonance) and provide the first explicitly reported manifestation of ghost resonance in 
a geosystem, i.e. of a mechanism which could be relevant for other systems with thresh- 
olds and with multiple states of operation. Our work enables us to explicitly simulate 
realistic probability measures of Dansgaard-Oeschger events (e.g. waiting time distri- 
butions, which are a prerequisite for statistical analyses on the regularity of the events 
by means of Monte-Carlo simulations). We thus think that our study is an important 
advance in order to develop more adequate methods to test the statistical significance 
and the origin of the proposed glacial 1470-year climate cycle. 
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1 Introduction 



Time series of North Atlantic atmospheric / sea surface temperatures during the last ice age reveal 
25 the existence of repeated large-scale warming events, the so-called Dansgaard-Oeschger (DO) events 
(Dansgaard et al. 1982| Grootes et al. 1993). In climate records from the North Atlantic region the 
events have a characteristic saw-tooth shape (Fig. [TJ: They typically start with a warming by up to 
10-15 K ( |Severinghaus and Brook| [T999"! |Leu enberge r~et al.||1999] l over only a few years/decades. 
Temperatures remain high for centuries/millennia until they drop back to pre-events values over a 
30 century or so. A prominent feature of DO events is their millennial time scale: During Marine 
Isotope Stages (MIS) 2 and 3, successive events in the GISP2 ice core were reported to be often 
spaced by about 1470 years or multiples thereof ( |Alley et al] |2001a| |Schulz| |2002[ |Rahmstorf| 
2003 1, compare Fig. [T] A leading spectral peak corresponding to the 1470-year period was found 



(Groo tes" and Stuiver| |1997| Yiou et al. 1997| l, and this spectral component was reported to be 
35 significant at least over a certain time interval (Schulz 2002). We note, however, that the statistical 
significance of this pattern is still under debate (Ditlevsen et al. 2007), in particular because of the 
lack of adequate nonlinear analysis methods. 

It was proposed that DO events represent rapid transitions between two fundamentally differ- 
ent modes of the thermohaline ocean circulation (THC) (Oeschge r et al.| |1984| |Broecker et al. 



40 |1985[ ), most likely corresponding to different modes of deep-water formation ( |Alley and Clark| 
|1999[ [Gano polski and Rahmstorf 2001). The origin of these transitions is also under debate: In 



principle they could have been caused by factors from outside of the Earth system ( |Keeling and 



|Whorf| [20001 [Riall [20041 |Clemens| [2003j |Braun et aL] [2003] l, but they could also represent internal 



oscillations of the climate system (Broecker et al. 1990 Sakai and Peltier 1997 van Kreveld et al. 



45 2000). Several nonlinear resonance mechanisms have been suggested in order to explain the char- 



acteristic timing of DO events, including coherence resonance ( Ganopolski and Rahmstorf 2002 



Timmerma nn et al.| |2003| |Ditlevsen et al.| |2005) and stochastic resonance (Alle y et aL"l |2001a|b"l 
Ganopolski and Rahmstorf 2002 Rahmstorf and Alley 2002). 



2 Spectrum of models 

50 DO events have already been simulated by a variety of models, ranging from simple conceptual ones 
to Earth system models of intermediate complexity (EMICs). Conceptual models are most suitable 
to perform large numbers of long-term investigations because they require very little computational 
cost. However, they are often based on ad-hoc assumptions and only consider processes in a highly 
simplified way. In addition to that, the number of adjustable parameters is typically large compared 

55 to the degrees of freedom in those models. This implies that seemingly good results can often be 
obtained merely by excessive tuning. Nevertheless, conceptual models can provide important help 
for the interpretation of complex climatic processes. 
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The gap between conceptual models and the most comprehensive general circulation models 
(GCMs), which are not yet applicable for millennial-scale simulations because of their large com- 
60 putational cost, is bridged by EMICs (Claussen et al. 2002| l. EMICs include most of the processes 



described in comprehensive models (in a more reduced form), and interactions between different 
components of the Earth system (atmosphere, hydrosphere, cryosphere, biosphere, etc.) are simu- 
lated. The number of degrees of freedom typically exceeds the number of adjustable parameters by 
orders of magnitude. Since many EMICs are fast enough for studies on the multi-millennial time 
65 scale, they are the most adequate tools for the simulation of DO events. 

3 The conceptual model 

The simple conceptual model which we use here is an extended version of the model described by 



Braun et al. ( 2005 | l (in the Supplementary Material of that publication). Here we use the model 



to demonstrate and analyse two apparently counterintuitive resonance phenomena (stochastic res- 
70 onance and ghost resonance) that can exist in a large class of highly nonlinear systems. Due to 
the complexity of many of those systems it is often impossible to precisely identify the reasons for 
the occurrence of these resonance phenomena. Our conceptual model, in contrast, has a very clear 
forcing-response relation as well as a very low computational cost and thus provides a powerful 
tool to explore these phenomena and to test their robustness. Furthermore, we describe and dis- 
75 cuss the applicability of the model for improved statistical analyses (i.e. Monte-Carlo simulations) 
on the regularity of DO events. In the following the key assumptions of the conceptual model are 
first discussed. In the Supplementary Material we then compare the model performance under a 
number of systematic forcing scenarios with the performance of a more comprehensive model (the 
EMIC CLIMB ER-2), compare Supplementary Information File and Supplementary Figs. 1-6. In the 
80 framework of the conceptual model we finally demonstrate and interpret two hypotheses that were 
previously suggested in order to explain the recurrence time of DO events, and we discuss how these 
hypotheses could be tested. 

3.1 Model description 

Our conceptual model is based on three key assumptions: 

85 1 . DO events represent repeated transitions between two different climate states, corresponding 
to warm and cold conditions in the North Atlantic region. 

2. These transitions are rapid compared to the characteristic life-time of the two climate states 
(i.e. in first order approximation they occur instantaneously) and take place each time a certain 
threshold is crossed. 

90 3. With every transition between the two states the threshold overshoots and afterwards ap- 
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proaches equilibrium following a millennial-scale relaxation process. This implies that the 
conditions for a switch between both states ameliorate with increasing duration of the cold 
and warm intervals. 

Our three assumptions are supported by paleoclimatic evidence and/or by simulations with a cli- 
95 mate model: 

1 . Since long, DO events have been regarded as repeated oscillations between two different cli- 
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mate states (Dansgaard et al. 1982 1. It has been suggested that these states are linked with 



different modes of operation of the THC (Oeschger et al. 1984 Broecker et al. 1985 ). This 
seminal interpretation has since then influenced numerous studies and is now generally ac- 
cepted (Rahmstorf 2002 ). Indirect data indicate that the glacial THC indeed switched between 
different modes of operation (Sarnthein et al. 1994 Alley and Clark 1999] > which, according 
to their occurrence during cold and warm intervals in the North Atlantic, were labelled stadial 
and interstadial modes. A third mode named Heinrich mode (because of its presence during 
the so-called Heinrich events) is not relevant here. 

2. High-resolution paleoclimatic data show that transitions from cold conditions in the North At- 
lantic region to warm ones often happened very quickly, i.e. on the decadal-scale or even faster 



( Taylor et al. 1997 Severinghaus and Brook 1999 1. The opposite transitions were slower, i.e. 
on the century-scale ( |Schulz| |2002 ), but nevertheless still faster that the characteristic life-time 
of the cold and warm intervals (which is on the centennial to multi-millennial time scale, com- 
pare Fig. 1). The abruptness of the shifts from cold conditions to warm ones has commonly 
been interpreted as evidence for the existence of a critical threshold in the climate system that 
needs to be crossed in order to trigger a shift between stadial and interstadial conditions (Alley 
et al. 2003| >. Such a threshold could be provided by the THC (more precisely, by the process 



of deep-water formation): When warm and salty surface water from lower latitudes cools on 
its way to the North Atlantic, its density increases. If the density increase is large enough (i.e. 
if the surface gets denser than the deeper ocean water), surface water starts to sink. Otherwise, 
surface water can freeze instead of sinking. The onset of deep-water formation can thus hinder 
sea-ice formation and facilitate sea-ice melting (due to the vertical heat transfer between the 
surface and the deeper ocean). A switch between two fundamentally different modes of deep- 
water formation can thus dramatically change sea ice cover and can cause large-scale climate 
shifts. Such nonlinear, threshold-like transitions between different modes of deep-water for- 



mation are at present considered as the most likely explanation for DO events ( Alley et al. 
[T999l[Ganopolski and Rahmstorf||200T] l. 

3. The time-evolution of Greenland temperature during the warm phase of DO events has a 
characteristic saw-tooth shape (Fig. fTJ. Highest temperatures typically occur during the first 
decades of the events. These temperature maxima are followed by a gradual cooling trend 
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over centuries/millennia, before the system returns to cold conditions at the end of the events. 
This asymmetry supports the idea that the system overshoots in some way during the abrupt 
warming at the beginning of the events and that the subsequent cooling trend represents a 
1 30 millennial-scale relaxation towards a new equilibrium ( |Schulz et al.| |2002| |Centurelli et al.| 

2006). We note that the time-evolution of Greenland temperature provides no clear evidence 



for an overshooting during the opposite transitions (i.e. from the warm state back to the cold 
one). This, however, is not necessarily in contradiction to our assumption: This lack of an 
overshooting in the temperature fields does not necessarily mean that the ocean-atmosphere 
135 system did not overshoot, since Greenland temperature evolution in the stadial state might 

have been dominated by factors other than the THC (respectively its stability), e.g. by Green- 
land ice accumulation, which would mask the signature of the THC in the ice core data. 

We will show later (in Sect. 3.3.) that the assumption of an overshooting in the stability of 
the system is in fact strengthened by the analysis of model results obtained with the coupled 

140 model CLIMBER-2. In that model the overshooting results from the dynamics of the tran- 

sitions between the two climate states: In the stadial state deep convection occurs south of 
the Greenland-Scotland ridge (i.e. at about 50 °N). In the interstadial state, however, deep 
convection takes place north of the ridge (i.e. at about 65 °N). The onset of deep convection 
north of the Greenland-Scotland ridge, which releases accumulated energy to the atmosphere 

1 45 (i.e. heat that is stored mainly in the deep ocean), in first place starts DO events in the model. 

This heat release leads to a reduction of sea ice, which in turn further enhances sea surface 
densities between 50 °N and 65 °N (e.g. by increased local evaporation and reduced sea ice 
transport into that area). As a result deep convection also starts between 50 °N and 65 °N, and 
much more heat can be released to the atmosphere. Without a further response of the THC the 

1 50 system would return quickly (within years or decades, i.e. with the convective time scale) to its 

original state. In CLIMBER-2, however, the changes in deep convection trigger a northward 
extension and also an intensification of the ocean circulation (i.e. an overshooting of the At- 



lantic meridional overturning circulation; compare Ganopolski and Rahmstorf (2001 1), which 
maintains the interstadial climate state since it is accompanied by an increase in the salinity 

1 55 and heat flux to the new deep convection area (at about 65 °N). In response to the overshooting 

of the overturning circulation, the system relaxes slowly (within about 1000 years, i.e. with 
the advective time scale) towards the stadial state. We note that the advective time scale cor- 
responds to the millennial relaxation time in our conceptual model. The model CLIMBER-2 
also supports the validity of our overshooting assumption during the opposite transition (from 

1 60 the warm state back to the cold one), as we will show in Sect. 3.3. We would like to stress that 

our interpretation of the processes during DO events is, of course, not necessarily true since 
we can only speculate that the underlying mechanism of the events is correctly captured by 
CLIMBER-2. 
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3.2 Model formulation 



165 We implement the above assumptions in the following way (compare Fig. |2j: First we define a 
discrete index s(t) that indicates the state of the system at time t (in years). Since we postulate the 
existence of two states, s can only take two values (s=l: warm state, s=0: cold state). We further 
define a threshold function T(t) that describes the stability of the system at time t (i.e. the stability 
of the current model state). 

1 70 Second we define rules for the time evolution of the threshold function T. When the system shifts 
its state, we assume a discontinuity in the threshold function: With the switch from the warm state 
to the cold one (at time t' in Fig. [2]) T takes the value Ao- Likewise, with the switch from the cold 
state into the warm one (at time t" in Fig. |2]i T takes the value Ai. As long as the system does not 
change its state the evolution of T is assumed to be given by a relaxation process: 
dT (T - B s ) 

175 — = -- (1) 

at t s 

(s labels the current model state, r s denotes the relaxation time in that state, B s is a state-dependent 
constant that labels the equilibrium value of T in each model state). These assumptions result in the 
following expression for the threshold function T: 

T(t) = (A s - B s ) ■ exp(-— —) + B s . (2) 

Tg 

180 Note that in the above expression the index s again denotes the current state of the model (i.e. s=0 
stands for the cold state and s=l for the warm one), S labels the time of the last switch from the 
warm state into the cold one, and 6i indicates the time of the last switch from the cold state into the 
warm one. 

Third we assume that transitions from one state to the other are triggered each time a given forcing 
1 85 function f(t) crosses the threshold function T. More precisely, we assume that when the system is in 
the cold state (s[t"] = 0) and the forcing is smaller than the threshold value (/[*" + 1] < T[t" + 1]) 
the system switches into the warm state (s[t" + 1] = 1). This shift marks the start of a DO event. 
Likewise, when the system is in the warm state (s\b'\ = 1) and the forcing is larger than the threshold 
value (f[t' + 1] > T[t' + 1]) the system switches into the cold state (s[t' + 1] = 0). That shift 
1 90 represents the termination of a DO event. If none of these conditions is fulfilled, the system remains 
in its present state (i.e. s[t + 1] = s[t]). 

To simplify the comparison of the model output with paleoclimatic records we further define a 
state variable S, which represents anomalies in Greenland temperature during DO events. For sim- 
plicity we assume that the state variable is equal to the threshold function: S(t) — T(t) (i.e. we 
1 95 assume that Greenland temperature evolution during DO events is closely related to the current state 
of the THC, respectively to its stability). We stress that this assumption is of course highly simpli- 
fied, because Greenland temperature is certainly not only influenced by the THC but also by other 
processes such as changes in ice accumulation during DO oscillation. However, this assumption is 



6 



not crucial for the dynamics of our model, since the timing of the switches between both model states 

200 is solely determined by the relation between the forcing function f and the threshold function T. This 
means that even if we included a more realistic relation between T and S, the timing of the simulated 
climate shifts would be unchanged and the model dynamics would thus essentially be invariant. 

Note that <5 and Si are not adjustable; they rather represent internal time markers. Thus, six 
adjustable parameters exist in our model as described here, namely Aq, Ai, Bo, B\, to and t\. Our 

205 choice for these parameters is shown in Table [T] With these parameter values the system is bistable 
(i.e. no transition is ever triggered in the absence of any forcing, since Bq < and Bi > 0) and 
almost symmetric. That means that the average duration of the simulated warm and cold intervals is 
almost equal. When compared with Greenland paleotemperature records this situation most likely 
corresponds to the time interval between about 27000 and 45000 years before present, during which 

210 the duration of the cold and warm intervals in DO oscillations was also comparable (Fig. 1). The 
model can, however, also represent an unstable (for B > and B x < 0) or a mono-stable system, 
in which the stable state is either the warm one (for Bo > and B\ > 0) or the cold one (for Bo < 
and B\ < 0); when compared to the ice core data this situation is closer to the time interval between 
15000 and 27000 years before present, since during that time the system was preferably in its cold 

215 state and the forcing apparently crossed the threshold only infrequently and during short periods of 
time. 

3.3 Comparison with a coupled climate model 

In order to test our conceptual model we compare its performance under a number of systematic 
forcing scenarios with the performance of the far more comprehensive model CLIMBER-2 (a short 
220 description of that model is given in the Appendix; a detailed description exists in the publication of 
IPetoukhovet al.| ( [2000l ». 

Analogous to Braun et al. ( 2005 ), we investigate the response of both models to a forcing that 



consists of two century-scale sinusoidal cycles. In the conceptual model, the forcing is implemented 
as the forcing function f. In the EMIC, the forcing is added to the surface freshwater flux in the 
225 latitudinal belt 50-70 °N, following panopolski and Rahmstorf| ( |2"00"T] l and |Braun etaij ( [2"0"0"5] l. This 



anomaly changes the vertical density gradient in the ocean and can thus trigger DO events. Switches 
from the cold state into the warm one are excited by sufficiently large (order of magnitude: a few 
centimetre per year in the surface freshwater flux into the relevant area of the North Atlantic) neg- 
ative freshwater anomalies (i.e. by positive surface density anomalies that are strong enough to 
230 trigger buoyancy [deep] convection), and the opposite switches are triggered by sufficiently large 
positive freshwater anomalies (i.e. by negative surface density anomalies that are strong enough to 
stop buoyancy [deep] convection). This justifies our choice for the logical relations that govern the 
dynamics of the transitions in the conceptual model (i.e. f(t) < T(t) as the condition for the switch 
from the cold state to the warm one, f(t) > T(t) for the opposite switch). 
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235 A detailed comparison between both model outputs is presented in the Supplementary Material. 
We here only summarise the main results: We find a general agreement between both models, which 
is robust when the forcing parameters are varied over some range (Supplementary Figs. 1-6). The 
conceptual model reproduces the existence of three different regimes (cold, warm, oscillatory) in 
the output of the EMIC and also their approximate position in the forcing parameter-space. By 

240 construction only the nonlinear component in the response of the EMIC to the forcing is reproduced 
by the conceptual model (this component represents the saw-tooth shape of DO events). A second, 
more linear component is not included in the conceptual model (this component represents small- 
amplitude temperature anomalies which are superimposed on the saw-tooth shaped events in the 
EMIC). In particular, the conceptual model very well reproduces the timing of the onset of DO 

245 events in the EMIC. The fact that our conceptual model, despite its simplicity, agrees in so many 
aspects with the much more detailed model CLIMBER-2 suggests that it indeed captures the key 
features in the dynamics of DO events in that model. 

We would like to stress that the output of the EMIC indeed supports our assumption of an over- 
shooting in the stability of the system during the transitions between both climate states: When 

250 driven by a periodic forcing (with a period of 1470 years), the EMIC can show periodic oscillations 
during which it remains in either of its states for more than one forcing period (i.e. for considerably 
more than 1470 years, compare Supplementary Fig. 4a). This implies that (at least in the EMIC) 
the conditions for a return to the opposite state indeed ameliorate with increased duration of the 
cold or warm intervals. If the thresholds in the model were constant (or gradually increasing with 

255 increasing duration of the simulated cold / warm intervals), in contrast, the duration of the cold and 
warm intervals during the simulated oscillation could never be longer than 1470 years: If a periodic 
forcing does not cross a constant (or gradually increasing) threshold within its first period, it never 
crosses the threshold, due to the periodicity of the forcing. 

4 Nonlinear resonance mechanisms in the model 

260 Strongly nonlinear systems can show complex and apparently counterintuitive resonance phenomena 
that cannot occur in simple linear systems. In this section we use our conceptual model to demon- 
strate and to discuss two of these phenomena, i.e. stochastic resonance (SR) and ghost resonance 
(GR). Since the explanation of the 1470-year cycle (and in fact even its significance) is still an open 
question, we further discuss how future tests could distinguish between the proposed mechanisms. 

265 4. 1 Stochastic resonance (SR) 

In linear systems that are driven by a periodic input, the existence of noise generally reduces the 
regularity of the output (e.g. the coherence between the input and the output). This is not necessarily 
the case in nonlinear systems: Excitable or bistable nonlinear systems with a threshold and with 
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noise, which are driven by a sinusoidal sub-threshold input, can show maximum coherence between 
270 the input and the output for an intermediate noise level, for which the leading output frequency of 
the system is close to the input frequency. This phenomenon is called stochastic resonance (SR) 
( |Benzi et al.| [19 82 ; Gammait oni et al.||1998| l. SR has been suggested to explain the characteristic 



timing of DO events (Alley et al. 2001a I, i.e. the apparent tendency of the events to occur with a 
spacing of about 1470 years or integer multiples thereof. It has further been demonstrated that DO 

275 events in the model CLIMB ER-2 can be subject to SR (Ganopolski and Rahmstorf 2002 ). 

Here we apply our conceptual model to reproduce these results and to reanalyse the underlying 
mechanism. We use an input that is composed of: (i) a sinusoidal signal with a period of 1470 years, 
(ii) additional white noise. Figures [3] and [4] show that for a suitable noise level the model can indeed 
show DO events with a preferred spacing of about 1470 years or integer multiples thereof. The 

280 reason for this pattern in the output is easily understandable in the context of the model dynamics: 
DO events in the model are triggered by pronounced minima of the total input (total input = periodic 
signal plus noise). These minima generally cluster around the minima of the sinusoidal signal, and 
the start of the simulated events thus has a tendency to coincide with minima of the sinusoidal signal 
(Fig. [3^). Some minima of the sinusoidal signal, however, are not able to trigger an event, because 

285 the magnitude of the noise around these minima is too small so that the threshold function is not 
reached by the total input. Consequently, a cycle is sometimes missed, and the spacing of successive 
events can change from about 1470 years to multiples thereof. 

Unlike the model CLIMBER-2 (which has a complex relationship between the input and the out- 
put and also a large computational cost) our conceptual model can be used for a detailed investigation 

290 of the SR, e.g. because the dynamics of the model is simple and precisely known and because proba- 
bility measures (such as waiting time distributions) can be explicitly computed. In fact, the resonant 
pattern in the conceptual model (Fig. [3} is due to two time-scale matching conditions: The noise 
level is such that the average waiting time between successive noise-induced transitions is compa- 
rable to half of the period of the periodic forcing, and also comparable to the relaxation times To 

295 and n of the threshold function (compare Fig. [£[5). This situation is different from the usual SR, in 
which thresholds (or potentials) are constant in time (apart from the influence of the periodic input). 



In the usual SR, only one time-scale matching condition exists (Gammaitoni et al. 1998 1, namely 
the one that the average waiting time between successive noise-induced transitions (i.e. the inverse 
of the so-called Kramers rate) is comparable to half of the period of the periodic forcing. 

300 In order to investigate the implications of the second condition we simulate histograms for four 
different scenarios in the conceptual model (Fig.|4|: 1. noise-only input, constant threshold (Fig.|4^); 
2. noise-only input, overshooting threshold (Fig. [4}}); 3. noise plus periodic input, constant threshold 
(Fig. [4]:); 4. noise plus periodic input, overshooting threshold (Fig. |4jl). We note that 3. corresponds 
to the usual SR, while 4. describes our overshooting stochastic resonance. As can be seen from 

305 the histograms, the existence of the millennial-scale relaxation process leads to a synchronisation in 



9 



the sense that the waiting times between successive events are confined within a much smaller time 
interval (about 1000-4500 years with the overshooting, compared to about 0-10000 years without 
overshooting). 

This confinement is plausible since the transition probability between both model states strongly 

310 depends on the magnitude of the threshold, which declines with increasing duration of the cold 
or warm intervals: When the standard deviation of the noise level is chosen such that the average 
waiting time between successive noise-induced transitions is comparable to the relaxation times To 
and t\, as in Fig.|4j the overshooting relaxation strongly reduces the transition probability for waiting 
times much smaller than the relaxation time (since the corresponding values of the threshold function 

315 are large) and increases the transition probability for waiting times of the order of the relaxation time 
or larger (since the corresponding values of the threshold function are considerably smaller). The 
probability to find an only century-scale spacing between successive events is thus small, because 
the corresponding transition probabilities are small. On the other hand, the probability to find a 
multi-millennial spacing is also small, because the states are already depopulated before (i.e. the 

320 probability to obtain lifetimes considerably larger than the relaxation time is almost zero). This 
explains why the possible values for the spacing between successive DO events are restricted to a 
much smaller range than in the usual SR (i.e. in the case with constant thresholds). 

This synchronisation effect is indeed not unique to the conceptual model: The output of the cou- 
pled model CLIMBER-2 shows a similar pattern (with possible waiting times between successive 

325 DO events of e.g. about 1500-5000 years or about 1000-3000 years, depending on the noise level; 
compare Fig. 4a-d in the publication of Ganopolski and Rahmstorf ( 2002| l). This similarity is of 



course not surprising, since the conceptual model is apparently able to mimic the events in the EMIC 
and since an overshooting in the stability of both states clearly also exists in CLIMBER-2 (compare 
Sect. 3.3). We note that in the GISP2 ice core data, DO events in the time interval 27000-45000 
330 years before present (which, as discussed in Sect. 3.2, is the best analogue to the "background cli- 
mate state" in our conceptional model, since the duration of the cold and warm intervals in the data 
is comparable in that interval) have spacings of about 1000-3000 years (compare Fig. 1) and were 
reported to cluster around values of either about 1470 years or about 2940 years (Schulz 2002). 



Because the SR mechanism could explain such a pattern (compare Fig. [4jl) it has originally been 
335 proposed. However, this mechanism requires a sinusoidal input with a period of about 1470 years, 
which has so far not been detected. 

4.2 Ghost resonance (GR) 

In linear systems which are driven by a periodic input, the frequencies of the output are always iden- 
tical to the input frequencies. This is not necessarily the case in nonlinear systems. For example, 
340 nonlinear excitable (or bistable) systems that are driven by an input with frequencies corresponding 
to harmonics of a fundamental frequency (which itself is not present in the input) can show a reso- 
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nance at the fundamental frequency, i.e. at a frequency with zero input power. This phenomenon, 



which was first described in order to explain the pitch of complex sounds (Chialvo et al. 2002 



Chialvo 2003 1 and later observed experimentally e.g. in laser systems (Buldu et al. 2003 i, is called 



345 ghost resonance (GR). GR and SR can indeed occur in the same class of systems, e.g. in bistable 
or excitable systems with thresholds. However, unlike SR, GR requires a periodic driver with more 
than one frequency. Although many geophysical systems might be subject to GR (since the relevant 
processes often have thresholds), the occurrence of this mechanism has so far not expressly been 
demonstrated in geoscience. 

350 Here we discuss a hypothesis that was recently proposed to explain the 1/(1470 years) leading 



frequency of DO events (Braun et al. 2005 i. The underlying mechanism of the hypothesis is in fact 
the first reported manifestation of GR in a geophysical model system. According to that hypothesis 
the 1/(1470 years) frequency of DO events could represent the nonlinear climate response to forcing 
cycles with frequencies close to harmonics of 1/(1470 years). Our conceptual model illustrates the 

355 plausibility of this mechanism: We use a bi-sinusoidal input with frequencies of 7/(1470 years) and 
17/(1470 years), i.e. with frequencies corresponding to the 7th and the 17th harmonic of a 1/(1470 
years) fundamental frequency, and with equal amplitudes. A spectral component corresponding to 
the fundamental frequency is not explicitly present in the input. Since the two sinusoidal cycles 
correspond to harmonics of the missing fundamental, the input signal repeats with a period of 1470 

360 years. For an appropriate range of input amplitudes, the output of the conceptual models shows 
periodic DO events with a period of 1470 years (Fig.|5]l. Unlike the input, the model output exhibits 
a pronounced frequency of 1/(1470 years), corresponding to the leading frequency of DO events 
and to the fundamental frequency that is absent in the input. This apparent paradox is explained by 
the fact that the two driving cycles enter in phase every 1470 years, thus creating pronounced peaks 

365 spaced by that period. Because the magnitude of these peaks results from constructive interference 
of the two driving cycles, it is indeed robust that a threshold process can be much more sensitive to 



these peaks than to the two original driving cycles (Chialvo 2003 1. 

The main strength of the GR mechanism is that - unlike the SR mechanism - it can relate the 
leading frequency of DO events to a main driver of natural (non-anthropogenic) climate variability, 
370 since proxies of solar activity suggest the existence of solar cycles with periods close to 1470/7 
(=210) years (De Vries or Suess cycle) and 1470/17 («86.5) years (Gleissberg cycle) (Stui ver and 



|Braziunas} 1993; Wagne r"er al.| |200"T] |Peristykh and Damon| |2003| ). So far, however, no empirical 
evidence for this mechanism has been found (Muscheler and Beer| |2006| l, nor has it been shown 
yet that changes in solar activity over the solar cycles are sufficiently strong to actually trigger DO 
375 events. 

In order to investigate the stability of this mechanism we further add a stochastic component (i.e. 
white noise) to the forcing. In this case the events are - of course - not strictly periodic anymore. 
Similar to the SR case, an optimal (i.e. intermediate) noise level exists for which the waiting time 
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distribution of the simulated events exhibits a maximum at a value of 1470 years, corresponding to 
380 the period of the fundamental frequency of the two input cycles (Fig.|6|l. In contrast to the SR case, 
in which a fairly simple waiting time distribution with a few broad maxima of century scale width is 
obtained (compare Fig. |4[l), we now find a much more complex pattern with a large number of very 
sharp lines of only decadal scale width. Since the waiting time distributions of both mechanisms 
are considerably different, it could - at least in principle - be possible to distinguish between both 
385 mechanisms by analysing the distribution of the observed DO events. In practise, however, this 
approach is complicated by the fact that only about ten events appear to be sufficiently well dated 



for this kind of analysis, and even their spacing has an uncertainty of about 50 years (Ditlevsen 



|et al. 2007 1, which is already of the same order as the width of the peaks in Fig. ^jp. We note that 



the mechanism that is described in Fig. [6] is known as ghost stochastic resonance (GSR), and its 
390 occurrence and robustness has already been reported before in other systems with thresholds and 



multiple states of operation (Chialvo 2003). At least in our system, however, this mechanism is 
even more complex than the other two types of resonance (SR, GR). It is beyond the scope of this 
paper to describe the GSR mechanism in more detail. 

4.3 Testing the proposed mechanisms 

395 The most direct way to test which of the proposed mechanisms - if any - provides the correct 
explanation for the timing of DO events would be to reconstruct decadal-scale density anomalies in 
the North Atlantic in connection with the events. This is not possible since even the most highly 
resolved oceanic records do not allow to reconstruct the variability of the glacial ocean on that time 
scale. Thus, only indirect tests can be performed. The identification of the postulated 1/(1470 years) 

400 forcing frequency, which has so far not been detected, would certainly give further support for the 
SR mechanism. And in order to support the GR mechanism, it remains crucial to demonstrate that 
century-scale solar irradiance variations are indeed of sufficiently large amplitude to trigger repeated 
transitions (with a preferred time scale of about 1470 years) between the two glacial climate states. 
This could be tested with climate models. 

405 An elegant and simple test is to make use of the observation that DO events in the Earth system 
model CLIMBER-2 represent the nonlinear response to the forcing, and that an additional (and much 
smaller) linear response is superimposed on the events. In the absence of any threshold crossing (e.g. 
in the Holocene, during which DO events did not occur) the response to the forcing, in contrast, does 
not show a strong nonlinear component. This suggests that Holocene climate archives from the North 

410 Atlantic region might be able to reveal what triggered glacial DO events. This approach has two 
major advantages: First, more reliable (e.g. better resolved and dated) records are available to solve 
this issue. Second, linear analysis methods can be used for that purpose, e.g. linear correlations. In 
the context of the GR mechanism, for example, the existence of a pronounced correlation between 
Holocene climate indices from the North Atlantic and solar activity proxies (reconstructed e.g. from 
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41 5 14 C variations in precisely dated tree rings) would be expected. Up to now at least one study exists 
that supports this prediction of a linear relationship between century-scale solar forcing and North 
Atlantic climate variability throughout the Holocene: Proxies of drift ice anomalies in the North 
Atlantic show a persistent correlation and a statistically significant coherency with "rapid (100- to 



200-year), conspicuously large-amplitude variations" in solar activity proxies (Bond et al. 2001 1, in 

420 accordance with the proposed GR mechanism. 

The most challenging test, however, is the direct analysis of the glacial climate records. We are 
convinced that one of the main difficulties in this approach is the high degree of nonlinearity of 
the events, which - according to our interpretation - has so far not been adequately addressed in 
many previous studies. For example, several attempts have already been made in order to investi- 

425 gate the 1470-year cycle by means of linear spectral analysis methods, and significance levels have 
commonly been calculated by assuming a red noise background. To us this assumption seems to 
be oversimplified, since the system responds at a preferred time scale even when driven by white 
noise (compare Fig. |4|)). We thus suspect that the significance levels obtained by this method are 
unrealistically high. We further think that the reported lack of a clear phase relation between solar 

430 activity proxies and DO events (Musc heier and Beer] |2006| l cannot rule out the idea that solar forc- 
ing synchronised DO events, since in the case of an additional stochastic forcing component (i.e. in 
the GSR case) the events are triggered by the combined effect of solar forcing and noise. Thus, the 
observed lack could also imply that only some of the events were in first place triggered by the Sun, 
whereas others were caused mainly by random variability (e.g. by noise). 

435 A new and promising approach, which is based on a Monte-Carlo method, has recently been pro- 
posed in order to test the significance of the glacial 1470-year climate cycle: Ditlevsen et al. ( 2QQ7\ 
define a certain measure in order to distinguish between different hypotheses for the timing of DO 
events, and they explicitly calculate the value of this measure for the series of events observed in 
the ice core. They then compare the calculated value with the values obtained by several hypothetic 

440 processes, e.g. by a random process (for which assumptions concerning the probability distribution 
of the recurrence times of the events have to be made). Significance levels are obtained from the (nu- 
merically estimated) probability distributions of the measure as generated by the considered process. 
Although we do not share their conclusions (because we think that more adequate measures can be 
chosen, which give considerably different results) we think that this approach is elegant because 

445 significance levels are not calculated based on linear theories. The method is thus also applicable to 
highly nonlinear time series. A major hurdle in this method is that for each considered process the 
probability distribution of the waiting times - which is unknown for almost all processes - somehow 
has to be specified. For example, Ditlevsen et al. use a simple mathematical (i.e. an exponential) 
distribution in order to mimic random DO events. In order to improve their novel approach, some 

450 method is thus needed to calculate waiting time distributions in response to any possible input. Com- 
prehensive models are not applicable, due to their large computational cost. Our conceptual model, 
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in contrast, is well designed for that purpose because it is combines the ability to mimic the complex 
nonlinearity of DO events as described by an accepted Earth system model with the extremely low 
computational cost of a very simple (zeroth order) model. We thus think that our work is an impor- 
455 tant step in order to develop improved statistical analysis methods which are able to cope with the 
extreme nonlinearity of DO events. 

5 Discussion and conclusions 

We here discussed DO events in the framework of a very simple conceptual model that is based on 
three key assumptions, namely (i) the existence of two different climate states, (ii) a threshold pro- 

460 cess and (iii) an overshooting in the stability of the system at the start and the end of the events, which 
is followed by a millennial-scale relaxation. These assumptions are in accordance with paleoclimatic 
records and/or with simulations performed with CLIMBER-2, a more complex Earth system model. 
In a couple of systematic tests we showed (in the Supplementary Material) that despite its simplic- 
ity, our model very well reproduces DO events as simulated with CLIMBER-2, whose dynamics 

465 is based on a (albeit reduced) description of the underlying hydro-/thermodynamic processes. The 
correspondence between both models thus strengthens our interpretation that the conceptual model 
can successfully mimic key features of DO events, and that these can be regarded as a new type of 
non-equilibrium oscillation (i.e. as an overshooting relaxation oscillation) between two states of a 
nonlinear system with a threshold. 

470 Although we discussed our model dynamics in the context of the (thermohaline) ocean circulation, 
our model does not explicitly assume that DO events are linked with changes in the ocean circulation: 
Threshold behaviour and multiple states exist in many compartments of the climate system (not only 
in the ocean, but e.g. also in the atmosphere and in the cryosphere). Our model thus cannot rule out 
a leading role of non-oceanic processes in DO oscillations. The millennial time scale of the events 

475 (which is represented in our model by the assumption of a millennial-scale relaxation), however, 
corresponds to the characteristic time scale of the thermohaline circulation and thus points to a key 
role of the ocean in DO oscillations. 

The main strength of our model is its simplicity: Due to the obvious relationship between forcing 
and response, the model can demonstrate why even a simple bistable (or excitable) system with a 

480 threshold can respond in a complex way to a simple forcing, which consists of only one or two 
sinusoidal inputs and noise. We applied our model to discuss two highly nonlinear and apparently 
counterintuitive resonance mechanisms, namely stochastic resonance and ghost resonance. In doing 
so we reported a new form of stochastic resonance (i.e. an overshooting stochastic resonance), 
in which the overshooting of the system leads to a further synchronisation effect compared to the 

485 usual stochastic resonance. Our study provides the first explicitly reported manifestation of ghost 
resonance in a geophysical (model) system. Since threshold behaviour and multiple equilibria are 
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not unique to DO events but exist in many geophysical systems, we would indeed expect that ghost 
resonances could be inherent in many geosystems and not just in our model. 

In addition to its applicability to demonstrate and interpret nonlinear resonance mechanisms, and 

490 to test their stability, we further illustrated the ability of our conceptual model to simulate probability 
measures (e.g. waiting time distributions, which are required in order to test the significance and 
the cause of the proposed glacial 1470-year climate cycle by means of Monte-Carlo simulations). 
Because it combines the ability to reproduce essential aspects of DO events with the extremely low 
computational cost of a conceptual model (which is up to 10 7 times lower than in the Earth system 

495 model CLIMBER-2), we think that our model represents an important advance in order to develop 
adequate nonlinear methods for improved statistical analyses on DO events. 



Appendix: Description of CLIMBER-2 

The Earth system model CLIMBER-2, which we used for our analysis, has dynamic components of 
the atmosphere, of the oceans (including sea ice) and the vegetation. Dynamic ice sheets were not 
500 included in our study. CLIMBER-2 is a global model with coarse resolution: For the atmosphere 
and the continents the spatial resolution is 10° in latitude, and 7 sectors are considered in longitude. 
The ocean is zonally averaged with a latitudinal resolution of 2.5° for the three large ocean basins. 
A detailed description of the model is given in the publication of Petoukhov et al. (2000). 



DO events in the model represent abrupt switches between two different climate states (stadial 
505 [i.e. cold] and interstadial [i.e. warm]), corresponding to two different modes of the THC: In the 
interstadial mode, North Atlantic deep water (NADW) forms at about 65 °N and much of the North 
Atlantic is ice-free. In the stadial mode, NADW forms at about 50 °N and a considerably larger area 
of the North Atlantic is ice-covered. We note that for the climatic background conditions of the Last 
Glacial Maximum (LGM) only the stadial mode is stable in the model whereas the interstadial mode 
510 is excitable but unstable (Ganopolski and Rahmstorf 2.00 1 1 > . Moreover, the stability of both modes 



depends on the actual climate state (e.g. on the configuration of the Laurentide ice sheet and on the 
freshwater input into the North Atlantic), and the stability properties of the system change when the 
background conditions are modified (more precisely, the system can be bistable or mono-stable). 
Transitions between both modes can be triggered by anomalies in the density field of the North 
515 Atlantic, for example by variations in the surface freshwater flux (since the density of ocean water 
increases with increasing salinity). In our study we thus implement the forcing as a perturbation in 
the freshwater flux (in the latitudinal belt 50-70 °N): We start the model with the climatic background 
conditions of the Last Glacial Maximum (LGM). Following earlier simulations ( |Braun et~aL 2005} 



we then add a small constant offset of 17 mSv (1 mSv = 10 3 m 3 /s) to the freshwater flux. For 
520 this climate state (which we label perturbed LGM) the THC is in fact bistable and DO events can be 
triggered more easily than for LGM conditions. This perturbed LGM state gives us the background 
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conditions for the model simulations as presented in this paper. 
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Tables 



Table 1. Parameters of the conceptual model. All parameters have the same values as in the publication 
of |Braun et aT] ( [2005) . Ao, Ax, B and Si are given in freshwater units (i.e. in mSv = milli-Sverdrup; 
1 mSv = 10 3 m 3 /s), since the conceptual model was originally designed to mimic the response of the THC to 
an anomaly in the surface freshwater flux. 



Parameter 


Chosen value 


Aq 


-27 mSv 


Ai 


27 mSv 


B„ 


-9.7 mSv 


Bi 


11.2 mSv 


T 


1200 years 


Tx 


800 years 
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Figure Captions 



Fig. 1. DO events as seen in the GISP2 ice-core from Greenland (Grootes et al. 1993, Grootes and Stuiver 
1997). The figure shows Greenland temperature changes over the interval between 10000 and about 40000 
years before present. DO events (0-10) manifest themselves as saw-tooth shaped warm intervals. Dashed lines 
are spaced by 1470 years. 

Fig. 2. Dynamics of our conceptual model. Shown is the time evolution of the model, in response to a forcing 
that is large enough to trigger switches between both model states. Top: Forcing f (black) and threshold function 
T (red). Bottom: Model state s (grey; s=0 corresponds to the cold state, s=l to the warm one) and state variable 
S (green). At time t" the forcing falls below the threshold function and a shift from the cold state into the warm 
one is triggered. With this transition, the threshold function switches to a non-equilibrium value (representing 
an overshooting of the system) and afterwards approaches equilibrium following a millennial-scale relaxation. 
At time t' the forcing exceeds the threshold function, and a transition from the warm state back into the cold 
one is triggered. With this transition, the threshold function switches to another non-equilibrium value and 
approaches equilibrium following another millennial-scale relaxation, until the forcing again falls below the 
threshold function and the next switch into the warm state is triggered. Note that the state variable S is chosen 
to be identical to the threshold function T. For convenience, discontinuities in T and S are eliminated by linear 
interpolation. T, S and f are normalised in the figure. 

Fig. 3. Stochastic resonance. The input consists of: 1. a sub-threshold sinusoidal signal with a period of 1470 
years and an amplitude of 4.5 mSv (about 40 percent of the threshold value Bi above which DO events occur 
in the model), 2. a random Gaussian-distributed signal with white noise power signature (standard deviation 
cr = 8 mSv) and a cutoff frequency of 1/(50 years). The cutoff is used since no damping exists in the model 
and it thus shows an unrealistically large sensitivity to high-frequency (i.e. decadal-scale or faster) forcing. A: 
Total input (black), periodic input component (grey), model output (green). Dashed lines are spaced by 1470 
years. B: Relative frequency to obtain a spacing of 1470 years ±10% (triangles) respectively ±20% (squares) 
between successive events, as a function of the noise level a. 

Fig. 4. Distribution of the spacing AT between successive events. The input in a and b consists of noise only, 
with a standard deviation a of 8 mSv (as in Fig.[3^). In c and d, a sinusoidal input component (amplitude = 4.5 
mSv, period = 1470 years; compare Fig. [3} is added to the noise. In a and c the threshold function is constant in 
each state (20 mSv in the warm state and -20 mSv in the cold one), while the overshooting relaxation assumption 
is used in b and d (with threshold parameters as shown in Table 1). Thus, c corresponds to the "usual" SR while 
d shows our "overshooting" SR. 

Fig. 5. Ghost resonance. Top: Forcing (black) and model response (green). Middle: Amplitude spectrum 
of the forcing. Bottom: Amplitude spectrum of the model response. We use two sinusoidal forcing cycles, 
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with frequencies of 7/(1470 years) and 17/(1470 years), respectively, and with equal amplitudes. These two 
cycles coincide every 1470 years and create peaks of particularly pronounced magnitude, spaced by exactly that 
period. Thus, despite the fact that there is no spectral power at the corresponding frequency (see middle panel), 
the forcing repeatedly crosses the threshold at those intervals. Consequently, the response of the conceptual 
model (i.e. the time evolution of the state variable S) shows strictly repetitive DO events with a period of 1470 
years (as indicated by the dashed lines, which are spaced by 1470 years). Despite the lack of a 1470-year 
spectral component in the forcing, the output shows a very prominent peak at the corresponding frequency. 

Fig. 6. Ghost stochastic resonance. The input consists of: 1. two sinusoidal forcing cycles with frequencies 
of 7/(1470 years) and 17/(1470 years), respectively, and with an amplitude of 8 mSv, 2. a random Gaussian- 
distributed signal with white noise power signature and a cutoff frequency of 1/(50 years), as in Fig. [3] In a, 
the relative frequency to obtain a spacing of 1470 years ±1% (triangles) respectively ±2% (squares) between 
successive events is shown as a function of the noise level a. B shows the distribution of the spacing AT 
between successive events (standard deviation of the noise: 5.5 mSv). 
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Fig. 1. DO events as seen in the GISP2 ice-core from Greenland (Grootes et al. 1993, Grootes and Stuiver 
1997). The figure shows Greenland temperature changes over the interval between 10000 and about 40000 
years before present. DO events (0-10) manifest themselves as saw-tooth shaped warm intervals. Dashed lines 
are spaced by 1470 years. 
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Fig. 2. Dynamics of our conceptual model. Shown is the time evolution of the model, in response to a forcing 
that is large enough to trigger switches between both model states. Top: Forcing f (black) and threshold function 
T (red). Bottom: Model state s (grey; s=0 corresponds to the cold state, s=l to the warm one) and state variable 
S (green). At time t" the forcing falls below the threshold function and a shift from the cold state into the warm 
one is triggered. With this transition, the threshold function switches to a non-equilibrium value (representing 
an overshooting of the system) and afterwards approaches equilibrium following a millennial-scale relaxation. 
At time t' the forcing exceeds the threshold function, and a transition from the warm state back into the cold 
one is triggered. With this transition, the threshold function switches to another non-equilibrium value and 
approaches equilibrium following another millennial-scale relaxation, until the forcing again falls below the 
threshold function and the next switch into the warm state is triggered. Note that the state variable S is chosen 
to be identical to the threshold function T. For convenience, discontinuities in T and S are eliminated by linear 
interpolation. T, S and f are normalised in the figure. 
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Fig. 3. Stochastic resonance. The input consists of: 1. a sub-threshold sinusoidal signal with a period of 1470 
years and an amplitude of 4.5 mSv (about 40 percent of the threshold value Bi above which DO events occur 
in the model), 2. a random Gaussian-distributed signal with white noise power signature (standard deviation 
a = 8 mSv) and a cutoff frequency of 1/(50 years). The cutoff is used since no damping exists in the model 
and it thus shows an unrealistically large sensitivity to high-frequency (i.e. decadal-scale or faster) forcing. A: 
Total input (black), periodic input component (grey), model output (green). Dashed lines are spaced by 1470 
years. B: Relative frequency to obtain a spacing of 1470 years ±10% (triangles) respectively ±20% (squares) 
between successive events, as a function of the noise level a. 
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Fig. 4. Distribution of the spacing AT between successive events. The input in a and b consists of noise only, 
with a standard deviation a of 8 mSv (as in Fig.[3^). In c and d, a sinusoidal input component (amplitude = 4.5 
mSv, period = 1470 years; compare Fig. [3} is added to the noise. In a and c the threshold function is constant in 
each state (20 mSv in the warm state and -20 mSv in the cold one), while the overshooting relaxation assumption 
is used in b and d (with threshold parameters as shown in Table 1). Thus, c corresponds to the "usual" SR while 
d shows our "overshooting" SR. 
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Fig. 5. Ghost resonance. Top: Forcing (black) and model response (green). Middle: Amplitude spectrum 
of the forcing. Bottom: Amplitude spectrum of the model response. We use two sinusoidal forcing cycles, 
with frequencies of 7/(1470 years) and 17/(1470 years), respectively, and with equal amplitudes. These two 
cycles coincide every 1470 years and create peaks of particularly pronounced magnitude, spaced by exactly that 
period. Thus, despite the fact that there is no spectral power at the corresponding frequency (see middle panel), 
the forcing repeatedly crosses the threshold at those intervals. Consequently, the response of the conceptual 
model (i.e. the time evolution of the state variable S) shows strictly repetitive DO events with a period of 1470 
years (as indicated by the dashed lines, which are spaced by 1470 years). Despite the lack of a 1470-year 
spectral component in the forcing, the output shows a very prominent peak at the corresponding frequency. 
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Fig. 6. Ghost stochastic resonance. The input consists of: 1. two sinusoidal forcing cycles with frequencies 
of 7/(1470 years) and 17/(1470 years), respectively, and with an amplitude of 8 mSv, 2. a random Gaussian- 
distributed signal with white noise power signature and a cutoff frequency of 1/(50 years), as in Fig. [5] In a, 
the relative frequency to obtain a spacing of 1470 years ±1% (triangles) respectively ±2% (squares) between 
successive events is shown as a function of the noise level a. B shows the distribution of the spacing AT 
between successive events (standard deviation of the noise: 5.5 mSv). 
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